#!/usr/bin/perl -w
##This script substitute reads sequence and quality score in the bam file. also modify the cigar string 
## 1st end fastq to forward strand, 1st end bam:  the same sequence
## 1st end fastq to reverse strand, 1st end bam:  reverse complement the fastq sequence
## 2nd end fastq to forward strand, 2nd end bam:  the same sequence
## 2nd end fastq to reverse strand, 2nd end bam:  reverse complement the fastq sequence

##not completed yet!!!!

## author: Yaping Liu  lyping1986@gmail.com 
## time: 2013-1-15

#Usage:  perl substitute_seq_in_bam.pl [option] sample.bam sample.fastq type_of_end(1:1st end; 2:2nd_end) > output.sam 

my $MERGE_BAM = "perl /home/uec-00/yapingli/code/mytools/perl/mergeBamWithSameSM.pl"; 
my $SAMTOOLS = "/home/uec-00/shared/production/software/samtools/samtools";

use Getopt::Long;
my $second_end="";
GetOptions(
	"second_end" => \$second_end,
);

my $bam=shift(@ARGV) or die;
my @fastqs=@ARGV or die;


my $forward_bam = $bam;
if($bam =~ /\.sam$/){
	$forward_bam =~ s/\.sam$/.fwd.sam/;
	if($second_end ne ""){
		$cmd="$SAMTOOLS view -h -S -F 16 -f 128 $bam > $forward_bam\n";
	}
	else{
		$cmd="$SAMTOOLS view -h -S -F 16 -f 64 $bam > $forward_bam\n";
	}
}
else{
	$forward_bam =~ s/\.bam$/.fwd.sam/;
	if($second_end ne ""){
		$cmd="$SAMTOOLS view -h -F 16 -f 128 $bam > $forward_bam\n";
	}
	else{
		$cmd="$SAMTOOLS view -h -F 16 -f 64 $bam > $forward_bam\n";
	}
}


system($cmd);

my $reverse_bam = $bam;
if($bam =~ /\.sam$/){
	$reverse_bam =~ s/\.sam$/.rev.sam/;
	if($second_end ne ""){
		$cmd="$SAMTOOLS view -h -S -f 16 -f 128 $bam > $reverse_bam\n";
	}
	else{
		$cmd="$SAMTOOLS view -h -S -f 16 -f 64 $bam > $reverse_bam\n";
	}
}
else{
	$reverse_bam =~ s/\.bam$/.rev.sam/;
	if($second_end ne ""){
		$cmd="$SAMTOOLS view -h -f 16 -f 128 $bam > $reverse_bam\n";
	}
	else{
		$cmd="$SAMTOOLS view -h -f 16 -f 64 $bam > $reverse_bam\n";
	}
}


system($cmd);

foreach my $fastq(@fastqs){
	open(FH,"<$fastq") or die;
my %seqs;
my %quals;
my $linecount_within = 1;
my $linecount_global = 1;
my $key="";
#if(0){
	while(<FH>){
	chomp;
	my $line=$_;
		if ($linecount_within == 1)
        {
            # Double check the format
            if($line =~ /^\@(\S+)\#\w+\/\d/){
            	$key=$1;
            	#print "$key\n";
            }
            else{
            	print STDERR "Incorrect FASTQ file $fastq\nLine ${linecount_global}: $line\nMod4 lines should start with \@\n"

            }
  
        }
        elsif ($linecount_within == 2)
        {
            
            	$seqs{$key} = $line;

            
        }
        elsif ($linecount_within == 4)
        {
            
            	$quals{$key} = &qual64_to_qual33($line);

            
        	$linecount_within = 0;
        	$key="";
        }
	$linecount_within++;
	$linecount_global++;
}
close(FH);
#}


my $cmd="";

##forward strand


open(FH,$forward_bam) or die;

while(<FH>){
	chomp;
	my $line = $_;
	if($line=~/^\@/){
		print "$line\n";
		next;
	}
	my @in = split "\t";
	my $key = $in[0];
	if(exists $seqs{$key}){
		$in[5] = length($seqs{$in[0]})."M";
		$in[9]= $seqs{$key};
		$in[10]= $quals{$key};
		#print "$in[5]\t$in[9]\t$in[10]\tlength($seqs{$in[0]})\t$seqs{$key}\t$quals{$key}\n";
	}
	
	my $out = join "\t",@in;
	print "$out\n";
}

close(FH);


##reverse strand


open(FH,$reverse_bam) or die;

while(<FH>){
	chomp;
	my $line = $_;
	if($line=~/^\@/){
		#print "$line\n";
		next;
	}
	my @in = split "\t";
	my $key = $in[0];
	if(exists $seqs{$key}){
		$in[5] = length($seqs{$in[0]})."M";
		$in[9]= &reverse_complement($seqs{$key});
		$in[10]= reverse($quals{$key});
		#print "$in[5]\t$in[9]\t$in[10]\tlength($seqs{$in[0]})\t$seqs{$key}\t$quals{$key}\n";
	}
	
	my $out = join "\t",@in;
	print "$out\n";
}

close(FH);

##combine


}

`rm $forward_bam`;
`rm $reverse_bam`;

sub reverse_complement{
	my $dna = shift @_;

	# reverse the DNA sequence
     my $revcomp = reverse($dna);

	# complement the reversed DNA sequence
        $revcomp =~ tr/ACGTNacgtn/TGCANtgcan/;
        #print "$dna\n$revcomp\n\n";
        
        return $revcomp;
}


sub qual64_to_qual33{
	my $in = shift @_;
	my @quals = split "", $in;
	my @quals33;
	foreach my $qual(@quals){
		push(@quals33,chr(ord($qual)+33-64));
	}
	my $out = join "",@quals33;
	return $out;
}